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ABSTRACT 

A  method  is  developed,  using  the  quadratic  form  of  canonical  un- 
coupled state  variables  as  a  Lyapunov  function  to  determine  the  switch^ 
ing  logic  for  a  quasi-optimum  control  of  a  non-linear  dynamic  system0» 
Analytio  arguments  are  presented  to  support  the  method,  which  is  capable 
of  being  extended  to  any  order  system  and  any  number  of  controls „  sub- 
jeot  to  certain  praotical  limitations  noted  in  the  analysis 0 

A  computational  scheme  for  determining  the  canonical  state  variables 
and  control  functions  is  presented  and  a  topological  interpretation  is 
made  of  the  choice  of  variables  used  for  the  control  logic  calculation., 

The  controller  based  on  the  method  described  is  applied  to  a  non~ 
linear  second-order  system.  Phase  trajectories  for  both  the  uncontrolled 
and  controlled  systems  are  obtained  by  means  of  a  digital  computer  simula= 
tion*  Various  aspects  of  the  theoretical  limitations  of  the  method  pre- 
sented are  investigated  and  the  experimental  results  are  analyzed  with 
respect  to  the  theoretical  predictions,, 

#The  system  is  postulated  to  be  described  by  a  system  of  differential 
equations  of  known  forme 
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lo     Introduction., 

The   stability  of  dynamic   systems  whose  equations  of  motion  are  such 
that  their   solutions  are  difficult   or   impossible  to  obtain  in  closed  form 
may  successfully  be  analyzed  by  means   of  Lyapunov"s  Second  Methodo 

We  have  a  given  dynamic   system  described  by  the  set   of  differential 
equations  as 

x     m     gjx);     £(0)      =     0   0  (1-3) 

Then£  referring  to  Fig„  (l-l),  the  origin*,  according  to  LaSalle  and 

5 
Lefschetz  iss 

stable  whenever  for  each  R  ^A  there  is  an  r  ^R   such 
that  if  a  path  (a  motion)  g*  initiates  at  a  point  x°  of 
the  spherical  region  S(r)  then  it  remains  in  the  spherical 
region  S(R)  ever  after 5  that  is,  a  path  starting  in  S(r) 
never  reaches  the  boundary  sphere  H(R)  of  S(R)g 


asymptotically  stable  whenever  it  is  stable  and  in  addition, 
every  path  g"*"  starting  inside  some  S(R«)9  R~\  o.,  tends  to  tl 


sry  patn  g"*"  starting  msiae  some  t>(Rc)9  R^ 
origin  as  time  increases  indefinitely? 


the 


unstable  whenever  for  some  R  and  any  rs   no  matter  how  smalls 
there  is  always   in  the  spherical  region  S(r)  a  point  x 
such  that  the  path  g     through  x  reaches  the  boundary  sphere 
H(R)0 


<-Ki&J 


Figo  (l-l)   Geometric  Interpretation  of  the  Various  Types 
of  Stabilityo 


Now  suppose  we  have  a  given  dynamic  system  described  by  the  set  of 
differential  equations  written  in  vector  form  as 

x  =  £(x,u)  (1-2) 

or 

x  =  Fx  ♦  Du  +  g  (x,u)  (1-3) 

where  F  is  a  constant  square  matrix,  D  is  the  control  distribution 
matrix,  u  is  the  control  function  vector,  and  £  (x,u)  contains  terms  of 
a  second  and  higher  power  in  x  and  u.  The  linearized  form  of  the  equa= 
tion  is 

x  =  Fx  4     Du  (1-4) 

and  can  be  said  to  adequately  describe  the  motion  of  the  system  for 
small  values  of  x. 

According  to  Kalman  and  Bertram  ,  a  Lyapunov  function  may  be  de- 
fined  as  some  scalar  function  of  the  state  variables  V(x)  with  the  fol- 
lowing properties:  (a)  V(x)  }  0,  V(x)(  0  when  x  £  Xg,  where  Xg  is  the 
equilibrium  state,  and  (b)  V(x)  ■  V(x)  s  0  when  x  ■  Xg. 

We  wish  to  apply  Lyapunov' s  Second  Method  to  the  design  of  a  con- 
troller so  as  to  return  the  state  vector  (x)  to  a  position  of  equilib- 
rium (xg)  in  some  optimum  manner.  No  loss  of  generality  will  result 
if  we  choose  (x_g)  to  be  the  origin  of  the  state  space,  and  henceforth 
the  equilibrium  position  will  be  considered  as  such, 

Lyapunov  stated  that  if  a  Lyapunov  function  V(x)  exists  (being 
positive  definite)  in  some  neighborhood  of  the  origin,  and  if  -V(x) 
is  also  positive  definite  in  some  neighborhood  of  the  origin,  the 
system  is  asymptotically  stable  . 

In  general  the  method  employed  will  be  to  make  V(x)  as  negative  as 
possible  in  order  to  drive  the  system  to  the  origin  as  quickly  as 
possible. 


The  Second  Method  of  Lyapunov  leaves  open  to  the  individual  inves- 
tigator the  choice  of  a  particular  Lyapunov  function  to  be  used  in  the 
analysis  of  a  given  system.  For  a  given  system  thens  there  are  possibly 
an  unlimited  number  of  such  functions •  The  quadratic  form  of  the  canoni- 
cal state  variables  as  the  particular  Lyapunov  function  employed  in 
the  design  of  the  controller  was  governed  by  the  simplification  obtained 
in  the  control  logic. 

Analysis  and  design  for  a  controller  of  a  linear,  stationary  system 
using  the  quadratic  form  procedes  in  a  straightforward  manner,,  In  the 
case  of  a  non-linear  system  however,  the  limitation  that  the  linear  ap- 
proximation is  sufficiently  accurate  only  within  some  region  surrounding 
the  equilibrium  position  poses  a  rather  severe  restriction  on  the  extent 
of  the  state  space  for  which  a  given  control  can  be  useful* 

In  order  tc  overcome  this  limitation  to  some  degree  the  controller 
was  designed  in  such  a  manner  that  the  control  logic  was  obtained  by 
successive  linear  approximations  to  the  actual  system  at  points  on  the 
system  trajectory. 

For  best  performance  of  the  system  the  time  involved  in  the  calcula= 
tion  of  the  control  logic  should  be  infinitesimal  so  that  there  would  be 
no  appreciable  time  delay  in  the  application  of  the  proper  control  func- 
tion.  That  is,  the  system  state  vector  would  not  have  moved  from  the 
point (x, ),  where  the  state  variables  were  sampled,  to  a  new  position  on 
the  trajectory,  (xo)  where  the  control  function  was  actually  applied o 

In  actual  practice  a  finite  calculation  time  is  of  course  necessary 
so  that  the  state  vector  will  always  have  moved  from  the  point  of  sampling 
before  the  new  value  of  control  function  can  be  appliedo 

For  plants  with  long  time  constants,  such  as  might  exist  in  chemical 
or  metallurgical  processes  or  ship  stabilization  systems,  presently 

3 


available  digital  computers  have  computation  speeds  adequate  to  success-, 
fully  provide  on  line,  real  time  control. 

Application  of  this  method  to  plants  with  short  time  constants  will 
depend  upon  increased  speeds  of  the  computers  and  the  development  of  more 
refined  methods  of  calculation* 


2.     Theory, 

If  we  have  given  a  Lyapunov  function  V(x)     a     x Px  we  may  assert 
that     the  equilibrium  state  xe     ■     0  of  a  linear  dynamic   system 

x     =     F^c  (2-1) 

is  asymptotically  stable  if  and  only  if  given  any  symmetric  positive- 
definite  matrix  Q,  there  exists  a  symmetric,   positive  definite  matrix  P 
which  is  the  unique  solution  of  the  set  of  linear  equations  described  by 

pTp     4     PF     =     _Q  .  (2-2) 

A  simple  procedure  for  calculating  the  control  logic  attributed  to 
R.W.   Bass     procedes  as  follows:      Choose  Q^>  0  arbitrarily.     P>0  may  then 
be  calculated  from  (2-2).      ||x||2p  will  then  be  a  Lyapunov  function  for  the 
system  (2-1).     We  now  choose  u(t)   so  as  to  make  V  (x)   as  negative  as 
possible.      Let 

V(x)     a     xTPx    =     ||x||2P,     then 

V(x)     »     xTPx     4     xTPx  .  (2-3) 

We  wish  to  consider  the  effect  of  control  on  V.     Taking  the  transpose 
of  (1-4)  we  obtain 

£T     -     xTFT     t     uTDT  (2-4) 

and 

V(x)      s     xVpx     +     xTPFx     +     uTDTPx     +     xTPDu 
or 

V(x)     s     xT(FTP     +     PF)x     *      2uTDTPx   . 

T 
Since  P  is   symmetric,   DP     ■     PD.      Simplifying  this  gives 

V(x)     ■     -xTQx     4     2uTDTPx 

or  v(x)     =     -||X||2Q     ♦     2iiTDT^  •  (2-5) 

*  /  \  T  T 

To  make  V(x)  as  negative  as  possible  we  want  the  term  2u  D  Px  to  be 

as  negative  as  possible.  Therefore  the  control  logic  is  given  as 


Ui(t)  a  -ai  sgn  (HTTx)±  (2-6) 

where  the  a^  are  constants  denoting  the  maximum  allowable  magnitudes  of 
each  individual  control  effort  u^.   It  is  obvious  from  (2-6)  that  in 
order  to  make  the  term  2u  D  Px  as  negative  as  possible  the  individual 
controller  u. (t)  must  always  exert  a  maximum  effort ,  and  in  such  a 
direction  that  the  term 

-ai  sgn  (DTPx)i 
is,  in  fact,  negative.  Thus  the  control  logic  becomes  merely  logic  for 
detecting  the  appropriate  switching  points  of  the  "bang-bang"  controller 
Ui(t). 

In  the  case  where  there  is  only  one  controller  the  vector  u  may  be 
of  the  form  such  that  u   ■   (0,0,0  .  «  •  u)»  The  system  equation  is 


then  of  the  form 


x 


0 
0 


~Dn  "bn-l 


0  .  .  .  0 
1 


1 
.  -b. 


u 


(2-7) 


where  the  b^  are  the  coefficients  of  the  characteristic  equation  of  the 
unforced  system,  taken  in  reverse  order. 

If  we  make  the  equivalence  transformation  from  the  physical  state 
space  to  a  canonical  state  space  denoted  by  the  matrix  equation 

y_  -      Gx  (2-8) 

then  by  differentiating  we  get 

£  s  G±  (2-9) 

and  equation  (1-4)  becomes 


We  choose  G  such  that 

G"lFG     =   _A.  (2-10) 

v/here  .A.  is  a  diagonal  matrix  with  the  same  eigenvalues  as  F" . 

For  the  system  with  one  controller  such  as  described  by  (2-7)  we 
may  further  choose  the  particular  transformation  matrix  G  such  that 

GT^Du    =    £ 
where  E*     =     (a,  a,   a,    .    .    .a) •     This  has  the  effect   of  weighting  the 
value  of   all  the  state  variables  equally  in  determining  the  sign  of  u(t). 
If  the  P  of   (2-2)   is  chosen  to  be  the  identity  matrix  I  and  the   =j\= 
matrix  substituted  for  F,   equation  (2-2)  becomes 

jSji     ♦     I  A.      -     -Q.  (2=11) 

If  all  the  real  parts  of  the  eigenvalues  of  F  are  negative,   so  then  are 
all  those  of  _A_ ,   and  Q  is  positive  definite.     Also 

Hx)     =     £T(JLT     +    A)l*     2ET£.  (2-12) 

The  first  term  on  the  right  side  of  (2-12)  is  again  negative  and  the 

control  logic  is  given  by  the  second  ternu  For  a  single  controller, 

we  have 

n 
u(t)  -  -a  sgn  (  E  Y±)   •  (2-13) 

i 

In  other  words  the  switching  function  is  such  that  the  control  must 
change  sign  whenever  the  sum  of  the  n  canonical  state  variables  changes 
sign. 

Consider  the  given  form 

2.    "  -A-  L    +   "  D—  • 

#It  is  always  possible  to  perform  the  equivalence  transformation  (2-10) 
provided  F  has  no  repeated  eigenvalues.  A  system  described  by  an  nth  order 
differential  equation  having  repeated  roots  may  be  transformed  into  a  par= 
tially  uncoupled  set  of  1st  order  differential  equations.  It  will  be 
assumed  that  the  systems  considered  in  this  paper  have  no  multiple  roots. 


If  we  take  the  LaPlace  transform  we  get 

sY(s)  -  j(0)  =   -A-Y(s)  +  (H-DuCs) 
or 

Y(s)(sl  --A-)  =  ^(0)  +   GT^uCs)  .       (2-14) 

Since  from  (2-6)  the  control  variable  should  assume  only  values  of  +  a 

©r  -  a_    ,  we  stipulate  that  it  is  a  constant  during  a  given  control 

period  (between  switchings).  We  then  have  the  LaPlace  identity 

(G-lDu)(s)  =  _L  . 
s 

Solving  (2-14)  for  a  given  component  y^  we  get 

Yi(s)(s     -\±)     -     yi(0)     +    J-L 


or 


Yi(8)     =     7Ii^L    >     _J1__     .  (2-15) 

•      (s  -Xi)  s(s  -XA) 


1 
Then  the  time  domain  solution  of   (2-15)   is 


7±W     -     yi(0)eXit £i-     (1  -eXit) 

■*-i 


or 


\±t  a. 


yi(t)     =    |yi(0)     4     —L\  e*1       -         *         .    (2-16) 
»  <^i   J  X  ^ 

An  interesting  result  becomes  apparent  from  this  solution  as  t  in- 
creases without  bound.  For  the  case  where  \^   is  negative  (a  consequence 
of  Q  being  positive  definite)  we  see  that 

y^t)  _*.     *  ;  t  -*»  oo 


This  result  indicates  that  although  this  form  of  controller  is  useful  in 
making  v(x)  more  negative  in  an  asymptotically  stable  system,  it  will  not 
result  in  V(x)  -+■  0  as  t  -+■  OO  if  u  is  a  function  of  the  type  seen  in 
Fig.  (2-1). 


i 

Ht 

"-i 

■-^- 

-  "  ai 

U'A 


Fig.  (2-1)  Control  Function  U^  As  An  Ideal  Relay. 

In  fact  it  will  result  in  a  condition  of  chatter-motion  around  the  origin 
as  the  controller  becomes  dominant  at  points  in  the  state  space  where  the 
yi<||e||:  where  ||e||  assumes  some  particular  small  value. 

This  objectionable  feature  may  be  overcome  in  the  case  of  asymptoti- 
cally stable  systems  by  choosing  the  control  function  of  the  form  shown 
in  Fig.  (2-2)  below  so  that  the  controller  is  effectively  disengaged 
within  some  region  arbitrarily  close  to  the  origin,  allowing  the  system 
to  reach  the  equilibrium  condition  in  an  unforced  manner. 


J 
♦*r 1 

H*i 

1 a> 

m 


Fig.  (2-2)  Control  Function  V^   As  A  Relay  With  Dead  Zone. 


If  the  uncontrolled  system  were  such  that  one  (or  more)  roots  con- 
tained positive  real  parts,  (in  this  case  -Q  can  no  longer  be  positive 
definite)  then  for  some  particular  points  in  state  space,  it  is  still 


possible  that 

v(x)    =    -|x2||q  . 

However  V(x)>0  and  remains   so» 

Consider  a  plant  governed  by  the  equations 

dx/dt     =     Fx     +     Du(t) 
where  the  control  variables  are  subject  to  the  constraints, 

|ui(t)|    <    aj  <   oo      (i     =     1,    .    .    .,  m)    •  (2=17) 

This  is  essentially  the  relay  or  saturating  servo  problem,. 
The  problem  is  to  return  every  initial  state  to  the  origin<> 

It  is  well  known  that   if  F  has  eigenvalues  with  posi- 
tive real  parts  then  there  are  some  states  which  cannot  be 
returned  to  the  origin  by  any  control  subject  to  the  con= 
straints   (2-17)4. 

If  a  controller  of  the  form  described   is  applied  to  the  unstable 

system,    it  is  apparent  from  (2-5)  that  control  may  be  maintained   so   long 

as  the   inequality  (2-18)   below  exists .^     That   is 

-||x||2Q     +     2uTDTPx      <     Oo  (2=18) 

From  the  solutions   (2-16)  the  behavior  of  each  y^(t)  for  0<t<ks 

where  k  denotes  the  switching  period  is  as  in  Fig.    (2=3)« 


Fig.   (2-3)     Plot  of  y^t)  vs.  Time  For  \x±     =     -ai  sgn  (EfyJ^ 
Stable  Eigenvalue. 
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->-  i 


a; 
3T; 


Fig.  (2-4)  Plot  of  yi(t)  vs.  Time  For  u^  * 
Unstable  Eigenvalue,  y^(0)>|  ai 

I  Xi 


.aj_  sgn  (E  v^, 


Fig.  (2-5)  Plot  of  JTi(t)  vs.  Time  For  u±     =  -ai  sgn  (E  x^i* 
Unstable  Eigenvalue,  yj(0)  <l  aj  I  . 

I  Xi  I 

In  the  event  that  the  system  distribution  matrix  D  were  such  that 

the  controls  were  uncoupled,  i.e.  for  each  y^  there  were  a  corresponding 

u^,  then  it  would  be  possible  to  design  a  controller  whose  individual 

elements  changed  sign  at  the  zero  crossings  of  the  individual  state 

variables. 
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From  Figo    (2-5)   it  can  be  seen  that  the  element  would   switch  at 
t     s     To      It   is  possible  to  calculate  T^     ■     ^(y^O))  from  (2-16)  for 
the  case  where    "^^>09     y^(0)>0,  u^  <0  a 

yi(t)    *   jyi(o)    -    ..,.a.1.Ae    i      ♦  — !^L  ,  let  yi(t)    s    0 


then 


{ 


yi(o) 


a*     ^     XtT  a 


^     I    e     x       « L 


X 


giving 


X.T  a. 


e  ■ 


i     I  AiyiCO)-*! 


and 


I 


Xi      ^  i  Vi<°>-»i 


(2-19) 


From  (2-19)    it   is  clear  that  for  y^O)     ■     09  T     s     _JL_    ln(l)     -     0S 
indicating  that  this   system  also  will  result   in  chatter-motion  around 
the  origino 

Seleotion  of  a  control  function  u^  of  the  form  of  Fig,   (2-2)  would 
be  unsatisfactory  in  this   case  because  the  uncontrolled   system  is  un- 
stable at  the  origin*     A  better  form  of  function  would  be  that  of  the 
form  of  a  saturating  amplifier  with  a  high  gain.      (See  FigQ    (2-6) )0 


(iW 


* 


Fige  (2-6)  Control  Function  uj  With  the  Characteristics  of  a 
Saturating  Amplifier* 
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Employing  LaSalle  and  Lefschetz5s  definitions  of  stability  des=> 
cribed  in  Section  19,  a  system  whose  eigenvalues  all  have  negative  real 
parts  is  asymptotically  stable  with  no  control s   merely  stable  with  a 
"bang-bang"  controller  (Fig.  (2-1)),  and  again  asymptotioally  stable  if 
the  controller  contains  a  dead  zone  as  in  Fig„  (2-2) 0   Given  an  arbit- 
rary H(R),  a  system  with  eigenvalues  lying  in  the  right  half  plane  is 
defined  as  unstable,  and  the  same  system  with  a  controller  of  the  form 
described  may  be  defined  as  stable  if  X  is  restricted  to  lie  within 
some  S(r)  whose  radius  r  is  of  some  necessarily  small  value.  The 
restrictions  on  the  allowable  values  of  r  are  necessarily  related  to 
the  constraints  on  the  magnitudes  of  the  control  functions  (2-17), 

Given  a  non-linear  system  described  by  the  equation  z  *  F(x)x, 
the  linear  approximation  of  the  non-linear  system  referred  to  by  Kalman 
and  Bertram  may  be  obtained  by  a  Taylor  series  expansion  of  the  coef- 
ficient matrix  F(x)  at  the  origin.  The  linearized^  constant,  matrix  F 

<iF°  (x) 

has  the  elements     1 ,  .      (The  matrix  F  is  the  Jacobian  of   F(x) 

d  xi 

with  respect  to  x  at  x  s  x(0)  s  0)  • 

Suppose  a  system  described  by  the  equation 

x  =  F(x)x  +  Du 
is  to  be  controlled  in  an  environment  where  the  disturbances  are  of 
such  a  magnitude  as  to  carry  the  state  variable  beyond  the  region  of 
validity  of  the  linear  approximation  given  by  (1=4).  Referring  to 
Fig.  (2-7),  this  condition  corresponds  to  an  x°  lying  outside  of  S(r)0 
By  expanding  F(x)  about  the  point  x°  and  obtaining  the  Jacobian  of 
F(x)  with  respect  to  x  at  x  s  x°  we  may  obtain  a  linear  approxima= 
tion  at  x° 

x  =  F  x  ♦  Du  .  (2-20) 
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Figo  (2-7)  Geometric  Interpretation  of  Successive  Linear 
Approximation  Process 0 


Control  logic  based  on  this  approximation  will  be  valid  while  the 
trajectory  remains  within  the  spherical  region  S(p)c  If  S(p  )  is  the 
region  of  validity  of  the  linear  approximation  due  to  the  n***  sampling* 
then  by  an  iterative  process  such  that  the  trajectory  from  Xn  to  Xn 
always  remains  inside  S(pn)  we  may  obtain  an  optimum  control  policy 
based  on  the  particular  Lyapunov  function  employed  for  the  complete 
trajectory  for  x°  to  xe  ■  0o 
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3.  Application  of  Theory  to  Design  of  a  Single  Controller. 

It  was  shown  in  Section  2.  that  the  control  logic  resulting  from 
the  choice  of  the  quadratic  form  of  the  uncoupled  state  variables  as 
the  system  Lyapunov  function  is  of  the  form 

u(t)  »  -a  sgn  J  ^pi     Y±\ 
when  the  system  has  only  one  controller •  Thus,  the  control  problem  is 
reduced  to  determining  the  sign  of  the  sum  of  canonical  state  variables 
and  insuring  that  the  control,  u(t),  assumes  the  opposite  sign©  The 
steps  involved  in  this  process  consist  of  the  following! 

lo  Measuring  the  physical  state  variables  (x)o 

2*     Calculating  the  canonical  state  variables  (yO  from  the 
measured  (x). 

3.  Forming  ^  yi» 

4*     Applying  the  proper  control  u(t)     ■     -a  sgnVy^o 
As  was   shown  in  Section  2»,  the  variables   (yj  may  be  calculated  by 
the  equation 

y_    s     G~Xx  .  (3-1) 

-1 
The  problem  then  reduces  to  evaluating  G  , 

It  may  be  shown  that  a  matrix  G  such  that 

G^FG  =  J\_  (3=2) 

where_/\-is  a  diagonal  matrix,  may  be  formed  by  calculating  the  eigen- 
vectors of  the  matrix  F  and  constructing  an  array  in  which  each  eigen- 
vector is  a  particular  column  of  the  array*  Once  a  G  has  been  obtained 
it  is  merely  necessary  to  calculate  the  inverse,  G"%  in  order  to  be  able 
to  calculate  the  state  variables  (^[)« 

The  elements  of  the  diagonal  matrix  _/\=are  the  eigenvalues  of  F0 
If  the  system  is  oscillatory  the  elements  of  _A-  will  necessarily  be 
complex,.  Since  the  elements  of  the  coefficient  matrix  F  of  the  original 
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differential  equation  are  all  real,  consideration  of  (3-2)  reveals  tl 
for  an  oscillatory  system  there  must  be  complex  elements  in  G  and  G". 

When  G  is  complex,  care  must  be  taken  in  evaluating  G*'*   If  we 
let  G  a  (feG  +  j$G,  and  attempt  to  form  G"1  -  (dlG)^1  «.  j(Ag)"1 
difficulty  arises  immediately0  Since  oomplex  roots  arise  in  conjugate 
pairs,  there  will  be  two  eigenvectors  (oolumns)  of  (C&G)  due  to  the 
identical  real  parts  of  the  conjugate  pairs 0  These  vectors  will  be 
identical,  and  therefore  linearly  related  so  that  ((&G)  will  be  singulars 
Similar  reasoning  shows  that  (<»G)  will  also  be  singular,  so  that  it  is 
impossible  to  obtain  (&G)"1  and  (<$  G)*1, 

Let  us  form  a  new  matrix  H  in  the  following  manner 0  If  the  order 
of  G  is  n,  the  order  of  H  will  be  2n  so  that  for  every  complex  element 
of  G  there  will  be  four  elements  of  H  arranged  in  a  square  pattem<>  The 
two  elements  on  the  principal  diagonal  of  the  square  have  a  value  equal 
to  the  real  part  of  the  corresponding  complex  element  of  G„  The  ab- 
solute values  of  the  remaining  two  elements  are  equal  to  the  imaginary 
part  of  the  corresponding  complex  element  of  Gs   the  element  in  the  lower 
left  hand  corner,  taking  the  positive  sign,  while  that  in  the  upper  right 
hand  corner,  the  negative  signc  Where  the  corresponding  element  in  G 
is  real,  zeros  appear  in  the  locations  in  H  corresponding  to  the  missing 
imaginary  parts  © 

If  now  we  form  H^,  then  we  may  obtain  the  elements  of  G"^-  from 
the  proper  locations  of  H~l0 

However  since  we  are  interested  in  obtaining  the  (yj  corresponding 
to  the  measured  (x),  the  simplest  procedure  is  to  form  the  2  x  n  sup- 
plementary  vector  (x)  which  has  alternately  the  real  and  imaginary 
parts  of  the  physical  state  variables0  Since  these  state  variables  may 
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be  measured  and   have  physical   significance,   they  may  only  be  real,    so 
that   (x)     consists  of  the  values  of  the  state  variables  alternating  with 
zeros.     Matrix  multiplication  of  E~l  on  the  right  by  (x)      results    in  a 
2  x  n  column  vector  whose   elements  are  the  real  and    imaginary  parts   of 
the   canonical  space  variables   (y_)  o 

EXAMPLE,     The  product  of  a  complex  2x2   square  matrix  and  a 
column  vector  with  real  elements 


yi 
y2 


R 


11  4  jIll     R12  +  JJ12 
R21  +  JJ21     ^2  4  ^22 


Then 


*1 

y2 


Rllxl     4     ^ll^l     *     R12x2     4     JI123C2 

R21xl     4     JI213C1     4     R22x2     4     JI22X2 
Rearranging  gives 

yi     ■     (R11X1     *     R12x2)     4     a(IH3Cl     4     I12x2) 

y2     a      (R21X1     4     hz^     4     ^I213C1     4     X22x2) 
Suppose   instead  we  form  a  4  x  4  matrix  of  the  form  H"    »     Then 


and 


<&7l 

Rn 

♦#11 

R12         -^12 

X^ 

Jtyi 

+jin 

Rll 

tfilZ     R12 

K 

0 

&y2 

R21 

*#21 

R22         '0I22 

x2 

iy2_ 

O^l 

% 

J*22       R22 

0 

(ky1    - 

R11X1 

4       hz*Z 

I 

A  yx    = 

J1!!^! 

4     PizH 

ay2     s     R21xl      4       R22x5 


&y2    ■    Ji 


2  in 


JI223C2 


Rearranging  gives 

yi  *  (Rllxl  +  R12x2>  +  O^l^l  +   x12x2) 


(3-3) 
(3-4) 
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y2  -   (R21X1  +  R22x2)  +  ^hlxl     +   X22X2^ 
which  is  identical  to  (3-3)  and  (3-4). 

It  is  evident  from  the  above  discussion  that  in  the  case  of  an  oseil° 
latory  system  the  transformation 

£  =  G^x  (3-5) 

may  map  a  given  physioal  state  variable  x^  from  a  point  on  the  real 
axis  to  some  point  in  the  complex  plane.  Thus  for  a  given  coordinate 
axis  in  the  physical  n- space  there  is  a  complex  plane  in  the  canonical 
space* 

Since,  however,  there  is  a  one  for  one  correspondence  between  points 
in  the  two  spaces,  a  control  which  succeeds  in  reducing  all  state  vari- 
ables to  zero  in  the  canonical  space  will  have  also  reduced  the  real 
state  variables  to  zero. 

The  computation  method  employed  to  determine  the  transformation 
matrix  G  was  based  on  a  computer  program  MA.TSUB,  programmed  by 
Louis  W.  Ehrlich  of  the  University  of  Texas,  which  calculates  eigen- 
values and  eigenvectors  of  a  matrix  up  to  order  50.  The  general  method 
employed  is  to  make  an  initial  guess  for  each  eigenvector  of  the  matrix 
F  and  to  converge  on  the  correct  value  by  an  iterative  scheme  due  to 
E.E.  Osborne  .  The  original  program  employed  a  guess  of  1.0  for  all  com- 
ponents of  the  eigenvectors ■   Since  it  was  assumed  that  the  system 
matrix  F  was  not  changing  at  more  than  a  moderate  rate,  a  modification 
was  made  to  the  program  so  that  the  most  recently  calculated  value  for 
each  eigenvector  was  selected  as  the  initial  guess  for  a  new  calculation 
thereby  speeding  the  computation.  The  output  of  the  MATSUB  program  was 
also  modified  so  that  the  eigenvectors  were  arranged  into  two  n  x  n 
matrices  corresponding  to  the  real  and  imaginary  parts  of  the  transfor- 
mation matrix  G, 
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Prom  the  two  n  x  n  matrices  the  2n  x  2n  matrix  H  was  formed  and  the 
inverse  taken.  The  canonical  state  variables  were  then  determined  and 
the  proper  control  calculated  and  applied » 

The  path  of  the  system  trajectory  in  the  physical  space  (x)  was 
then  calculated  for  a  selected  time  interval  between  sampling  instants 
by  the  Runge-Kutta-Gill  method.  Instantaneous  values  of  the  elements 
of  the  coefficient  matrix  F  »  F(x)  were  calculated  for  the  new  sam- 
pling point  and  the  procedure  repeated  until  the  trajectory  in  the 
canonical  space  was  within  a  pre-set  epsilon  region  of  the  origin  or 
a  given  time  had  elapsed. 

Initially  the  ideal  case  was  simulated  by  holding  the  solution  of 
the  system  trajectory  at  the  sampling  point  until  the  new  calculation 
of  the  control  function  was  made  and  the  control  applied 0  Simulation  of 
the  non-ideal  case,  a  finite  calculation  time,  was  affected  by  applying 
the  control  calculated  at  the  (n-l)   samplings  at  the  time  of  the  n^h 
sampling.  This  is  equivalent  to  a  controller  continually  applying  a 
new  control  as  soon  as  it  can  perform  the  sampling  and  control  cal- 
culations. 

Simplified  flow  charts  of  both  the  ideal  and  non-ideal  case  siraula=> 
tions  are  presented  in  Appendix  I.  The  computer  program  in  FORTRAN 
language  (for  the  ideal  case)  is  also  included  in  Appendix  I« 
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4.     Simulation  of  Control  of  a  Non-Linear  Second  Order  System. 
A.     Background. 
The  non-linear  VanderPol   equation 

X     -     (1  -x2)  x     +     x    -     0  (4-1) 

■was   selected  as   the    system  on  which  to  evaluate  the   controller.     Al- 
though it   is   only  a   second   order  equation  it  was  felt  that    since  it 
is   quite  non-linear    it  would  be  a  fair  test   of  the   controller  and  the 
two  dimension  phase   space  has  the  advantage   of  making   graphical  analysis 
more  simple. 

The  equation  exhibits  a  stable  limit   cycle  of  the  general  shape 
shown  in  Fig.    (4-1). 


^   X 


Fig.    (4-1)     Phase  Portrait  For  an  Uncontrolled  VanderPol  Equation. 
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Referring  to  (4-1),  when  x  ■  0,  the  system  is  most  unstable.   For 
this  point  the  roots  are  +.5  ±     j.866.   It  was  determined  from  trajec- 
tories of  the  uncontrolled  system  obtained  by  digital  computer  analysis 
that  when  the  system  was  in  the  limit  cycle  the  maximum  excursion  of 
displacement  was   ~   2.1925.   This  value  was  substantiated  by  an  analog 
computer  simulation.   From  (4-1)  the  eigenvalues,  or  roots,  of  the  system 
when  the  displacement  reaches  its  maximum  value  in  the  limit  cycle  were 
calculated  to  be  -3.5235  and  -.2835  so  that  the  locus  of  roots  of  equa- 
tion (4-1)  during  a  complete  circuit  of  its  stable  limit  cycle  is  shown 
in  Fig.  (4-3). 


.      \ 

-X~\J  J   . 

1 

i  +  j* 

.5 

+i 

.866 

-3.5235 

-.2835 

+  CT 

^_ 

-3 

-X 

-il 

.5 

-3 

.866 

Fig.  (4-3)  Root  Locus  of  the  Uncontrolled  VanderPol  Equation 
in  Stable  Limit  Cycle. 


Examination  of  Fig.  (4-1 )  or  equation  (4-1)  reveals  that  the  system  be- 
comes unstable  whenever  |x|  <  1.0  and  becomes  stable  again  whenever 
|x|>  1.0.  Thus  the  equation  presents  the  controller  with  the  problem 
of  controlling  a  system  which  contains  at  various  times,  stable  real 
roots,  stable  complex  roots,  and  unstable  complex  roots. 

It  is  informative  to  observe  that  a  stability  analysis  of  the 
VanderPol  equation  by  Lyaounov' s  Second  Method  yields  precisely  the 
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same  results  as  the  more  familiar  methods  of  analysis  such  as  the  root 
locus. 

Rewrite  (4-1)  in  vector  form.  Then 


or 


and 


i     *      [-?        (l-x*)]    £  (4-2) 

xx     =     Xg  (4-3) 

x2     =     ~xl     +     (l-xj)xg   .  (4-4) 


Choose 


then 


V(x)     =     ||x||     -     x*    4     4 


V(x)     -     2X2XX     +     2x2x2  (4-5) 

and  from  (4-3)  and   (4-4), 

v"(x)     =     2xlx2     +     2xg(-x^     +     (l-x^Xg) 
or 

V(x)  =  2x|(l-x*)  (4-6) 

so  that  V(x)  is  positive  definite  but  -V(x)  is  not  for  jx,|<l.  There- 
fore the  region  near  the  origin  (-l<x<l)  is  unstable.  This  is  pre- 
cisely the  result  obtained  by  observing  that  the  root  locus  passes  into 
the  right  half  plane  whenever  |x]J<l. 
B.  Procedure. 
It  was  decided  to  initially  investigate  the  action  of  the  control 
with  no  time  delay  between  calculation  of  the  control  and  its  application 
in  order  to  verify  that  the  method  of  control  and  the  simulation  pro- 
cedure (described  in  Section  3.)  were  feasible.   First  a  family  of  tra- 
jectories of  the  system  with  no  control  effort  was  obtained  from 
various  initial  conditions  chosen  inside  and  outside  the  limit  cycle 
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(Graphs  A-l  through  A-8).   Trajectories  were  then  obtained  using  the  same 
initial  conditions  with  the  control  system  activated  (Graphs  B-l  through 
B-8).  The  sampling  rate  was  set  at  .3  seconds  and  u  was  set  at  1*0.  An 
arbitrary  epsilon  region  surrounding  the  origin  was  established  so  that 
the  trajectory  was  considered  to  have  reached  the  origin  when  the  sum 
of  the  canonical  state  variables  squared  was  less  than  .001.  Whenever 
the  trajectory  entered  the  epsilon  region  of  the  origin  the  solution 
was  terminated. 

Effects  on  the  trajectories  due  to  varying  sampling  rates  (Graphs 
C-l  through  C-5)  were  then  investigated,  and  finally  two  trajectories  were 
obtained  for  the  system  simulating  the  effect  of  finite  calculation  time 
(Graphs  D-l  and  D-2).  These  trajectories  all  were  started  from  an 
initial  condition  close  to  the  uncontrolled  limit  cycle  trajectory. 

The  simulation  program  included  a  sub-roatine  for  the  graphical  pre- 
sentation of  the  system  trajectories.  The  graphs  noted  above  are  pre- 
sented in  Appendix  II. 
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D,  Discussion  of  Digital  Computer  Simulation* 

Comparison  of  the  trajectories  of  the  controlled  and  uncontrolled 
systems  starting  from  the  same  initial  conditions  (Graphs  A-l  through 
A-8  and  B-l  through  B-8)  clearly  indicates  that  the  control  system  with 
a  control  effort  =  1.0  and  a  sampling  rate  of  .30  seconds"  is  capable 
of  driving  the  system  into  the  unstable  region  surrounding  the  origin  and 
maintaining  it  within  some  region  of  lesser  extent  than  that  enclosed  by 
the  system's  uncontrolled  limit  cycle.  Referring  to  Fig.  (l-l),  if  the 
radius  of  S(R)  were  established  for  this  system  as  R  ■  .5,  then  the  con- 
trolled system  might  be  defined  as  stable  while  the  uncontrolled  system 
would  be  defined  as  unstable.   Graph  B-l,  initial  conditions  x  ■  .2, 
x  =  .2,  illustrates  the  system  trajectory  in  the  region  close  to  the 
origin.   The  chatter-motion  mode  of  operation  is  readily  apparent. 

The  theoretical  analysis  of  Section  2.,  predicting  that  by  employ- 
ment of  this  type  of  controller,  an  unstable  system  could  be  made  stable, 
but  not  asymptotically  stable,  and  that  a  chatter-motion  mode  of 
operation  would  result  around  the  origin,  is  substantiated. 

Examination  of  Graphs  C-l  through  C-5  reveals  that  an  increase  in  the 
control  system  sampling  rate  causes  the  trajectory,  having  once  arrived  in 
the  neighborhood  of  the.  origin,  to  remain  thereafter  within  a  smaller  re- 
gion of  the  origin  than  the  same  system  operating  with  a  low  sampling  rate. 

Curves  B-3  and  C-3,  both  trajectories  for  the  identical  initial  condi= 
tions  and  system  parameters,  demonstrate  remarkably  different  trajectories, 
each  of  which  eventually  arrives  in  a  neighborhood  close  to  the  origin. 

#The  time  interval  between  successive  calculations  of  the  control 
function.   The  Runge-Kutta  integration  of  the  trajectory  employed  an 
interval  of  0.03  seconds. 
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Examination  of  the   calculated   data  revealed  that    in  the  case  of 
Graph  D-3  the  initial  values   of  the  canonical   state  variables  were 

yi     =      .75     -     j.1561 

y2  =  .75  ♦  j.1561 
while  in  the  case  of  Graph  C-3  the  values  were 

yx  -   .75  -j.1561 

y2     ■      .75     -j.1561   . 
Since  the  control  function  u     ■     -1.0     sgn  (    5Z0&y^     +     ^2  A  Ya)  i 
it   is  evident  that  u  assumed  a  different   sign  in  each  case.     That  this 
was   so  can  be  determined  by  examination  of  the  diverging  paths  of  the 
trajectories  in  the  two  cases. 

It   seemed  evident  that  the  more  direct  trajectory  should  be  the 
correct   one,   however  the  possibility  was  raised  that  there  might   not 
be  a  unique   solution  for  the  control  function.     Six  additional  runs 
each  exactly  coincided  with  the  more  direct  trajectory  B-3,   indicating 
that  oerhaps  an  error  in  calculation  led  to  the   selection  of  the  op- 
posite value  of  control  function  rather  than  there  being  two  possible 
solutions  for  the  control  function. 

The  computation  routine  had  been  programmed   so  that  the  inverse 
transformation  matrix  was   printed   at  each  sampling  point.      It  was  dis- 
covered that  the  transformation  matrices   in  the  two  cases  were 
not  identical. 

INITIAL  INVERSE  TRANSFORMATION  MATRIX  FOR   GRAPH  B-3 

.35820E-10  .64051E-00  .50000E-00  .40032E-00 

-.64051E-00  .35820E-10  -.40032E-00  .50000E-00 

-.35820E-10  -.64051E-00  .50000E-00  -.40032E-00 

.64051E-00  -.35820E-10  .40032E-00  .50000E-00 
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INITIAL  INVERSE  TRANSFORMATION  MATRIX  FOR  GRAPH  C~3 


.37180E-08  .64051E-00  o50000E=00  C40032E=00 

-.64051E-00  .37180E-08  -o40032E~00  .50000E-00 

o50000E-00  .40032E-00  .65670E-09  .64051S-00 

-.40032E-00  e50000E-00  -.64051E-00  <>65484E-09 


In  order  to  determine  if  two  solutions  were  possible,  the  two  in- 
verse matrices  were  checked  against  the  original  matrix  to  see  if  both 
were  valid 0 

The  original  matrix  was  calculated  by  hand  in  the  same  manner  that 
the  computer  was  programmed  to  do0 

The  procedure  is  illustrated  below. 
The  initial  F  matrix  was  I  and  the  initial  eigenvalues  were 

-0625  *  j.7806.  Then  for  the  first  eigenvalue 

X, 


L?  -U26J  |_*y 


(-.625     +     j.7806) 


y 


or 


*2     =     (-.625     *     je7806)x1 
Xj     =     (-.625     -     j.7806)x2 


so  that 


x«     -     le0 

Xj  =  -o625  -  j.7806  „ 
In  a  similar  manner,  for  the  second  eigenvalue 


*2 


1.0 

-.625  ♦  j.7806  . 


Then  the  2n  x  2n  matrix  becomes 


0625 

-.7806 

-.625 

O7806 

e7806 

-  v625 

-.7806 

-.625 

1 

0 

1 

0 

0 

1 

0 

1 

When  this  matrix  is  multiplied  by  the  inverses,  B-3  and  C-3,  only  the 
inverse  B-3  results  in  the  identity  matrix,  demonstrating  that  the 
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calculation  resulting  in  C-3  is  not  an  independent  correct  solutions 

Comparison  of  the  false  trajectory  with  the  correct  one  indicated 
the  possibility  of  there  being  at  least  two  other  instances  of  erroneous 
calculations  besides  the  one  at  the  origiro  These  are  indicated  by  the 
abrupt  changes  in  direction  (outward)  of  the  trajectory  at  the  points 
x  =   ,65,  y  =  0.6,  and  x     s  -.5,  y  s   .45  .  Calculations  similar 
to  the  one  carried  out  above  verified  the  fact  that  erroneous  determina- 
tions of  the  correct  control  function  had  been  made* 

The  similarity  of  the  trajectories  of  curves  C-l„  C-2  and  D»l  in 
the  fourth  quadrant  raises  the  interesting  comparison  between  the  control 
obtained  by  Lyapunov's  method  and  the  method  of  optimum  control  employing 
the  calculation  of  the  optimum  switching  lines  in  negative  time  as  dis- 
cussed  by  I0  Flugge-Lotz  and  H0Ao  Titus  0 

We  may  speoulate  that  there  is  some  optimum  switching  line  for 
this  system,  emanating  from  the  origin  in  a  manner  similar  to  that 
in  Fig.  (4-4), 


*~X 


Fig.  (4-4)   System  Optimum  Switching  Line, 
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The  optimum  trajectory  for  this  system  would  then  proceed  from  the 
initial  condition  to  its  first  intersection  with  the  optimum  switching 
line,  change  the  sign  of  its  controller^,  and  proceed  along  the  optimum 
switching  line  to  the  origin. 

It  appears  evident  that  the  method  of  switching  logic  discussed  in 
this  paper  selects  a  switching  time  too  soon  in  the  cases  of  curves  C~l 
and  C-2,   In  the  case  of  curve  D-2  the  system  appeared  to  be  in  a  state 
of  indecision.  The  switchings  at  points  A  and  C  undoubtedly  would  have 
brought  the  trajectory  initially  closer  to  the  origin  than  the  one 
finally  selected  at  B.  The  switchings  at  B  and  D  served  only  to  hinder 
the  system  response.   Investigation  of  the  computer  data  again  indicated 
erroneous  calculations  made  for  switchings  at  points  Bj>  C  and  D,  although 
the  resulting  decision  at  C  proved  to  be  a  correct  one. 

Graphs  D-l  and  D-2,  showing  the  simulation  of  the  system  operating 
under  non-ideal  conditions,  again  substantiate  the  fact  that  higher 
sampling  rates  generally  produce  more  satisfactory  operation  of  the 
system.   It  is  interesting  to  note  that  the  controlled  system  operating 
with  a  time  delay  of  .3  seconds  eventually  settles  into  a  stable  limit 
cycle  of  approximately  one  half  the  size  of  the  uncontrolled  system. 

The  effect  of  the  controller  on  the  excursions  of  the  system  roots 
as  the  state  point  travels  along  a  trajectory  for  a  particular  initial 
condition  is  demonstrated  in  Figs.  (4-5)  and  (4-6) „  The  numbers  indicate 
time  in  seconds  to  reach  the  position  indicated,. 

The  root  locus  for  the  uncontrolled  system  starting  from  initial 
conditions  x  ■  .5,  x  =  .5  is  illustrated  in  Fig.  (4-5).  The  roots 
attain  their  most  stable  positions  after  the  trajectory  has  settled  into 
its  stable  limit  cycle. 
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12.7 
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7.3 


12.7 


Fig,  (4-5)  Root  Locus  for  Uncontrolled  System*  Initial 
Conditions  x  ■  05,  x  =  #5» 

As  the  controlled  trajectory  approaches  the  region  of  chatter- 
motion  surrounding  the  origin,  the  root  locations  depart  very  little 
from  their  most  unstable  position  (  +  ,,5  ±  j,866)  because  maximum 
instability  occurs  when  x  s  0,  The  most  stable  root  positions  on 
the  controlled  system  root  locus  occur  early  in  the  solution  while 
the  trajectory  is  relatively  far  from  the  origin. 


Fig.  (4-6)  Root  Locus  for  Controlled  System,  Initial 
Conditions  x  =   .5,  x  s   .5. 
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5,  Conclusions . 

The  operation  of  the  controller  in  regulating  a  simulated  non- 
linear system  described  by  the  VanderPol  equation  supports  the  theo- 
retical predictions  developed  in  the  design*  The  solution  obtained  by 
this  method  is  seen  to  be  less  than  optimum  because  the  logic  does  not 
cause  the  control  to  switch  at  the  optimum  switching  curve.   It  is 
apparent  from  the  figures  that  the  sampling  rate  for  determining  the 
eigenvalues,  and  hence  changes  in  control  action,  is  critical  for  this 
to  be  an  effective  control  design  procedure.  With  the  fine  sampling 
the  system  disturbances  were  rapidly  and  effectively  controlled* 

There  exist  today  very  few  direct  techniques  for  the  control  of 
non-linear  systems.  The  procedures  studied  here  provide  such  a  tech- 
nique.  The  method  has  the  advantage  of  being  applicable  to  non- 
linear time  varying  systems  of  any  order.   It  does  however  require  a 
digital  computer  in  the  control  link. 
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APPENDIX  I 

FLOW  CHARTS  AND  FORTRAN  LANGUAGE  PROGRAM 
FOR  DIGITAL  COMPUTER  SIMULATION 
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START 
(INITIALIZE) 

1 

t 

CALCULATE 
EIGENVECTORS 
OF  COEFFICIENT 
MATRIX  F 

* 

FORM 
INVERSE 
MATRIX 

* 

CALCULATE 

z 

^y^  is 

5     ^s.      YES 

STOP 

J[    NO 

CALCULATE  AND 
APPLY  CONTROL 

A 

CALCULATE  NEW 
TRAJECTORY  AND 
X  AT  END  OF 
SAMPLING  PERIOD 

OLD  EIGENVECTORS 

ARE   INITIAL 

GUESS  FOR  NEW 

EIGENVECTORS 

* 

FORM  NEW 
COEFFICIENT 

A 

MATRJ 

[X  F 

Fig.  1-1   FLOW  CHART  FOR  IDEAL  CASE  SIMULATION, 
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START 
(INITIALIZE) 

i 

\ 

CALCULATE 

EIGENVECTORS 

OF  COEFFICIENT 

MATRIX  F 

* 

FORM 
INVERSE 
MATRIX 

* 

CALCULATE 

/^    IS   N.        YES 

STOP 

<.  t,  ^s            *~ 

OLD  EIGENVECTORS 

ARE  INITIAL 

GUESS  FOR  NEW 

EIGENVECTORS 

J  so 

CALCULATE 
CONTROL 

t 

FORM  NEW 
COEFFICIENT 

MATRIX  F 

s'               ^\^      YES 

^    INITIAL     \__£> 

SAVE  CALCULATED 

CONTROL 

SET  CONTROL 

a       0.0 

I 

r 

\^loop  r/ 

T  NO 

CALCULATE  NEW 
TRAJECTORY  AND 
X  AT  END  OF 
SAMPLING  PERIOD 

SAVE  CALCULATED 

CONTROL 

SET  CONTROL  = 

PREVIOUS  CONTROL 

1 

\ 

Fig.  1-2   FLOW  CHART  FOR  FINITE  CALCULATION  TIME  SIMULATION 
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..JOB    GOLDS    '   ■       OPVAN    1 
PROG!;  n     ] 

DIMENSION      ..:  (  20,  20)  ,  AI  (">  0.20  )  .  OR  t  70    ">c\     ^Tfin    on\ 
lXRl20)tXDOT{20)fXl    900    fvlf 900), VALUES    |G|ft?0,20,, 
2DUB kC,HQ) ,DUBINV(io?UO * .YEXPfin)    'xfyp I M     ' 
3X2U00i,Y2?100NX3}l60)    Y3(100)?x1m00    ^y4(    6o) 
7nnn    £KH0^,ARMAA|'GRfGI, VALUES  ; 

C7000MREAD    51D.MMfiEfF;U,Ej  r.DELTA.TF, ALRS.ALIS.  IXRCJI, J-UHMI 

C        ,E       -       COEFFICIENT    GF    SECOND    TERM    PF    VAlsinF^POi     enu 

C         tHi>      =       RADIUS    OF    EP  FGION   HF    'ihtptm 

r         Ifoc=    ,,T?cE    LIlMiI    F0R    SOLUTION 

I      aggfi    ^SEIGEHSffi8i,gSL15LS?I«!*fft"855B*UMe   °F   ARCI™^   OFFSET    FRO, 
C      5"iJF9«Mf         '$§02 ??§&f.?5§!e    VARIABLES 
1201    ^ggMAT(  I5^ '■FIO-O/TFIO.O) 

'20VrS^{{f l&b  wV\u  5IH-JUU  WB*HJU  UWOUi./ 
ri2o$RF4RAtti&K2?if8?EiE?f  g"i?0E5V C0EFF  HATRIX 

1205  FORMAT(lHO)  '  c-"3' 

ill  IIftS*i?5i^&*g?t 

PRINT  1200 

LINH  =  i+ 

TO  =  0.0 

NUMPTS=0 

LAG  =  0 

DO  200  J  =  j-MM 

DO  200  I  =  1,MM 

GR(I,J)  =  UO 


|*    Gin,  j)  =o.p; 

200    AIU.J)  =o.p' 

1   SSV?IilAb?uiaTlSpE*8SW^8f5GftN8fiB1gfN  'EC0TRS  0F  AR  AN0  F0RMS 

nOOO  CALL  'gSUBCMM  .  ALRS  :ALl5|AND  GI 


DO    30    1=1, MM 
M    =    2*1 
K   *   M  -1 

30LXEXPN(M??^IiR5°RMS    D0UBLE    SI2E    TRANSFORMATION    MATRIX    CALLED    DUB 

DO    31     1=1, MM 

M    =    ?*I 

K    =    M-l  • 

DO    31    J    =    1,MM 

L=    2*J    -1 

N  =    2*  J 

DUB    IK,L)    =    GR(I,J) 

DUB    (M,NJ    =    GRU,J 
_.    .DU6(K,N)    =-GHI,     ) 
3.1    DUB    (M,L)    =    GI(     ,j) 
*   2 

CALL    GAUSS    3     ( NN    ,     .1E-06.,     DUB    ,    DJBINV         kfvi 

GO    T0F(1?2)!k£?RSE    °F    DLlg  DUBINV    INV    '    KEY) 

STOP 

?    f8SWI,U9H    GR    MATRIX    SINGULAR) 
I     L.UI\l  IN U c 

FORMS  CANONICAL  STATE  VARIABLES  YEXP  (CONTAINS  BOTH  REAL  AND  IMAG.  TERMS 


YEXP(I)=  0.0 


AO 


' 


DO       2    J    =  1,NN 
32    YJn£(I,n=,vYEXPm     +    DUBINVd.J)     *    X=XP(J,1) 


,R  =  CO 


DO  100  3  I  =  1,NN 
5  SQUARE  OF  CANONICAL  RADIUS  V£ 
o  YSQR  =  YSQR  v  Y        *  YEXP(I) 
66  -  LINE)  1400,  1401,  1 1+ 0  1 
1  '*  C  0    INT  1200 


F 
100 


CTCR 


]UC1r     \42C2'T'X^^^R(2),{iVALUES(I5J),I  =  l,2},J=l,2),(YEX^:}7I-; 
'1203 

|UDUBINVU,J),J«lfNN>fI-lfNN) 

1402  till!    1200^  !^2,U03flM3 
LINE  =  4 

GO  TO  1404 
1*403  PRINT  1205 


- 


•«-    3 


C       TESTS    FCR    CLOSENESS    TO    0RIG1 
IS8S    i^0°S   :   JSgRM005.1303.1303 

C 

1 


CALCULATION  OF  CONTROL  FUNCTION 
DO  10  06  I  =  1 
1006  YRSUM  =  YRSUM  +  YEXP(I) 
UOP   =  -  SIGNF  I, YRSUM) 

1SI  ^!^[pfT;     mWW'2))-]°-E-b>    "2.702.T01 

7o0  VALUES (2 , 1 }=  -VALUES ?2, 15 

70S  gSWfSii'2'-?  "VA,-UES<2.2) 

t  TDEL  =  T  +  DELTA 

J304  JFCTDEL  -  T  -  DT).  1300  ,  1301  ,  1301 

300  IF  TDEL  -  TF)  1302,777,777 

1302  AR<  MM,  1)    «j.t    F 

GO^S'fooBJ    *    {K°    '    XRCn     *      XR(1)) 

1303  CONTINUE 


,^,      '  21202  ,. 
1202       7j2^:-TR 


PRINT  ill 202 

41  202  FORMAT!  25H*s 
-    .   PRINT  1202 ', 

"  7/7  CONTINUE 


INSIDE  EPSILOM///) 
RADIUS  VECTOR///) 


,  YSQR 

C ALL  X GRAPH  f  NUMPTS , X 1 , Y 1 , 8 } 
GO  TO  7000 

13re^R^  FOR  SAMPLING  PERIOD 

1405  k\l\    hlhHE)     1U°5'  ,4°6' U06 


LIME  =  4 
1406  PRINT  1202,T,XRt 1 ) , XR { 2 ) 
LINE  =  LINE  +  1 
I  Ft  SCO  -  NUMPTS)  1304,750,750 
750  NUMPTS  =  I     rS  +  1    f3U» 'DU 

rsj  -■  xr<  i) 

Yl{   1PTS)  =  XR( 2) 

ro  n 

END 


SUBROUTI  (M,  j  ,  XR  7DT  ,  UOP) 

JI§8!^tfWl5&,xpoT(20,•  XCI20)'  C(U1 


1 )  =0.0 
C(2)  =C5 
Cf3)  =  0.5 
C(41  =1.0 
00  105  0   I  =1,4   • 
TC  =  T  ♦  C  ( I  )  *  DT 
DC  1051   J  =  KM 
!051  XC(J)  =  XR(J)  4-  C(I)  *  AK(  T-]  ,j) 
CALL  DERIV  iTC    ,  XC  ,  XDOT  J  UOP 


M  ) 


a 


1052 


+2.0*  AK(2,J)  ♦  2,C*AK(3,J)  +  AK( 4, J ) ) /6. 
XR(20>,  ARC  20^20)  ,AI(20,20),GR(20,20),GI(2G, 


*  XR(1)  +AR(M,2)  *  XR(25  +  UOP 


922 


519 


917 
811 


523 
9 


403 
C   NEX 
C,    P 

i     504 

4 

5 
10 


11 

12 

13 

106 

109 


1U 


1 

I  =  1 ,  N 

J  *  1 ,  N 

CR(I,J) 
CR( I, J) 
CI Clf J) 
CHI, J) 
TO  IA 
TO 


I 


TO  ID 
530  TO  IA 
40  TO  IB 
523  TO  IC 


1050   J  =1,M 
AK  (I.J)   =  OT  »  XOOT{ J) 

DO  105  2   J  =1  ,M 

J)  =XR(J)  +(AK(1,J) 
END 

SUB ROUT  I  ME  DERIV  (T 
DIMENSION  XDOT  (20)  , 
0)  ,VALU£S(2t 20) 

N  AR  , A  1,GR,GI, VALUES 
)T( 1)  =  XR12) 
XDOT( 2)  =  AR(M, 1 ) 

.ROUTINE    GSUB     (M    ,     ALRS     ,     ALIS     )  rt  „       nn    „MV 

DIM.  IN    AR(20-20)  tAU20,20)  rBR(20,  20)  ,BI ( 2 0*20) , CR( 20. 20 >  • 

CI (20,20),XR(20) ,X1 (20) ,YR{20),YI <20 ) , VALUES ( 2, 20 ) ,GR(20,20) , 
GI( 20,20) 

COMMON  CR  ,CI tGRiGIt VALUES 
N  =  M 
LOOP  = 
DO  5)9 
DO  519 
BRC It  J) 
ARC  It  J) 
Bid, J) 
AI  (  I,  J) 
ASSIGN  917 
ASSIGN  811 

-•  M 
GO  TO  555 
GO  TO  ID 
ASSIGN  914 
ASSIGN 
ASSIGN 
ASSIGN 
ISL=-1 

GO  TO  92   '  !■ 
ISL=0 
ALR  =  ALRS 
ALI  =  ALIS 
IT=1 

DO  504  1=1, M 

T  TWO  INSTRUCTIONS  DESIGNATE  INITIAL 
RECEDING  EIGENVECTOR  SOLUTION 
XR{ I )=GR{ I, LOOP) 
XI  (  I)=GI  (  W600?) 
DO  5  1*1. N  I 
AR( I, I )=ARU, I )-ALR 
AH  I,I)  =  AI(  I,  I  )-ALI 
I  J=l 
BIG=0. 
DO  13  I=1,N 
YRC I )=0. 
YI{ I)=0. 

on  ii  j=i?n 

YRC  I  i=YR(IH-AR(  I  , J)*XR (J )-AI U , J ) *XI  (J) 

YI  (  I)=YI  (  I  )+AI{  I  ,J)*XRU  )  +  AR{  I  ,J)*XI  {  J) 

AM=YR( I)*YR( I )+YI(I )*YI( I) 

IF  (AM-BIG)  13,13,12 

BIG=AM 

JJ=I 

CONTINUE 

IF (BIG)  109,106,109 

ICT^IOOO 

J  J  =  N 

ISL  =*  1 

GO  TO  92 

RQNR=0. 

RQNI=0. 

ROD=0. 

DO  14  1=1, N 

IR=RQNR+XR( I ) *YR{ I ) +XI ( I ) *YI  ( I ) 
RQNI=RQNI+XR( I)*YI ( I )-XI ( I )*YF  ( [  ) 
R3D  =  RQD  +  XR( I) *XR ( I ) +XI ( I ) »X I v I  - 
Ay,UR  =  RQNR/RQD 


GUESSES  FOR  EIGENVECTOR 


12 


/RQD 

IUR+AiMUI*AMUI 
81     TS=0. 

DO    15    I  =  1  ,  N 
150TS=TS+( YRJ  I  )-AMUR*XR(  I  )  +  AMUI*XI<  I  )  )**2+ 
*XHI)-AMUI*XR(I))**2 

%t     XR ( I  )  =  { YR ( J J ) *YR (I ) +Y I ( J  J  5  * Y I ( I ) ) / BI  G 

16  XI  m  =  {YR{JJ)*YI(I)-YI  .       R(I))/BIG 

A  K  i  J  J  }  —  j  o  0 

XI { JJ)=0.0 

IF  (TS  /  RQD  -  10.E-4)  20,20,18 
IF { IJ  -  20   19,20,20 
9  IJ=IJ+1 

;  rn  10 

20  ICT=IJ 

T2    =     IJ    +    20 
R=AMUR+ALR 
.I=AMUI+ALI 
I  —  N 
DO    310    I=1,N 
,       ARC  I, I)=AR(ItI)-AMUR 
310    AI  (  I,  I  )=AI  (  I  ,  D-AMUI 
GO    TO    29 
09    DO     ICO     I=1,N 

(1,1  )  =  AR(I,I)-ALR 
100    AI(  1,1 )=AI C I , I J-ALI 

2  9    U=IJ*1 
535    DO    27    1  =  2, 1^ 
1.1  =  1  —  1 
DO    27    J=I, IM1 

21  FM=AR{ I, J)*ARl I , J ) +A I ( I , J ) *AI { I , J ) 

R<  J  »  J  WAR  ( J  t  J )  +AI  (J,  J  )  *A  I  (  J,  J ) 
„       '  I)    24,24,22 

T1=AR(J,k] 
T2=AI(JtK) 
AR{  J,  K)=AR(|(,K) 
[  J  ,  K )  =  A I  {  f ,  K  J 
AR(I,K)=T1 

23  MU,K)=T2 
T1=XR(J) 
T2=XI(J) 

( J)=XRi I) 
XI( J)=XI (II 
XR(I)-T1 
XI ( I )=T2 
T  1  =  f  ;-s 

=  SM 
S  M  =  T 1 

24  Ir'    (SM)     25,27,25 

25  IF     (FM)    90,27,90 

90  8?"CfA.5iI.,H!*  J)+Aici,j)»Ai(j,jn/sM 

DO    ^o    K=J,MM 
„        ARn,KJ=AR(I  ,K)-RR«AR{  J,KWRI*AI  { 

26  AI(  I,K)=AI(I  ,K)-RR*AI  {  J  ,  K)  -R  I  *AR  { 
AR { I , J ) =0. 

- [ I, J)=0. 
( ! 5=XR( I )-RR«XR(J )+RI*XI( J) 

27  ^!N^,n-RR*XUJ,-RI*XR,J' 

GO  TO  IA 
530  SMALL=1000. 
DO  28  K=1,MM 
T  K  K  =  K 

I]=^{,5*K)*AR{K,K)  +  AnK,K)*AI(K,KJ 

_  _  ^     ir     ii\j      /r>0,r52»750 
?50    IF     (T1-SMALL)    751,23,28 
751    SMALL=T1 
K 


J  ,  K  ) 


23    CONTINUE 

GO    TO    ID 


13 


752 
30 


974 
753 


32 

33 


31 


41 


42 
43 


45 

46 

49 

47 


60 : 

1  16 

755 
50 


55 

51 

52 

G3 


53 
93 

501 
40  0 
402 


L  KK 

USD    753,30,30 
■=  1 

r=2000 

974    1=1, 
XRCD-0.0 

I)=0.0 
YR(  IZ)  =  1.0 
YI(IZ}=C.C 

12 
>1,C 
IF    (IZ-MM)    33,32,33 
111  =  2 

■ 
DO    31    I^IZZ  , ; 
YR{ I )=0o 
YI ( I)=0. 
IZZ=MM-I Z+2 
IF     (IZ-1)    95,49,95 

=  1 
BIG-O. 

I=IZZ,MM 
II=MM-I«-1 
^  1 1  -e- 1 

<;:;  =  o. 

SI=0. 

IF     ll-l)    42,44.4  2 

>:0    43    K=KK 

SR  =  SR+AR (  I  I , K } *YR ( K ) - A  I (  I  I  ,  K  )  * Y I  ( K ) 

SI=SI+AR(II.K)*YI (KJ+AI (II,K>*YR{K) 

T1=AR(  11,11  i*AR{  II,II)+AHII,II)*AI(  II,  II) 

! 1 1 )  =  ( AR( 1 1 , 1 1 ) » { XR { 1 1 } -SR ) +A I ( 1 1,  1 1 ) * ( XI ( II ) -SI ) ) /T 1 
YI(  II  )  =  (  ARf  II,ID*CXI(  ID-SI  )-AI  (II,  ID*(XR{  ID-SR)  )/Tl 
AM=YR (II) »YR (  1 1 ) + Y I ( I  I ) *  Y I( I  I) 
IF     (AM-8IG).    46,46,45 
J  J  =  I  T 
BIG=- 
CONTINUE 

4  7    1=1 , MM 

I  J  =  ( YR ( J J  5 *YR { I ) +YI { J  J ) *Y I ( I )  ) / BI  G 
XHI}  =  (YRUJ)*YI(I)-YICJJ)»YR(I))/BIG 

!  J  J  )  =  1  .  0  . 


I  JJ)=0»0  :.  - 
92  OG  6C1  1=1  » N 
DO  601  J  =  1?N 

:  I,  J)=BR{,£  ,J) 
AI(1,J)=BI(I,J> 
IF  (ISL)  755,50,60 
GO  TO  IC 
ALR=0, 
ALI=0. 

=  0.0 
DO  52  1=1, N 
YR( 1  )=0» 
YI( I )=0. 
GO  51  K=1,N 

YRC  I  )=YR( I )+AR( I,K)  *XR(K  )-AI ( I  ,K)*XI  (K) 
YIU)=YI(D+AR{I,K)ftXI(K)  +  AI(I,K)ftX^(K) 
ALR=ALR+XR( I )«YR( I)+XI!I)*YIII) 
ALI  =  ALI+XR(  D*YI  (I)-XI(I  )»YR(I  ) 
SUM=SUM+XR( I) »XR( I ) +XI ( I )»Xi ( I ) 
A!.";  =  alr/Su:-] 
=ALI/SUM 
AM=ALR»ALR+ALI*ALI 
TS  =  0* 

DO  5  3  1  =  1 p  M 

T1=YR( I )-ALR»XR( i)+ALI*XI( I) 
T2  =  YI  (  I  )-Al.R«XI  (  I  J-ALIftXR(I  ) 
TS=TS+T1*T1+T2«T2 

IF  (TS  /  SUM  -  10.E-14)  60,60,301 
IF  (IJ-MIT2)  99, 400, MOO 
IF  (IT  -  3)  402,990,402 

R  =  ALR  +  .1 


-:-.;' <- 


60 
63 


61 


AL 

IS 

IF 
V] 
72 
XR 
XI 


63 


62 
65 


XI 

DO 

T2 

AR 
AI 
AR 
AI 
DG 
7? 
12 
AR 
A I 
AR 
A  T 


'  =  AL 

70  h 
.  =  0 
S  {  1 
S(2 

I  JJ 
XI  (JJ 

J  J  I 
J  J  . 
N)=T1 
-T2 
6£  K  = 


iLR 


I  +  .1 


yLOOP)   = 

,LQOP)  =  ALI 
)  61,65,61 
) 
) 

UN) 


AR(  JJ 
=  AIUJ 

J  o  ?  K  / 
JJ,K) 
N  ,  K )  = 
K}  = 


66 


6C0 
700 


701 
702 


- 
DO 
DO 
AR 
Ai 
DO 
DO 

81 

DO 
DO 

A  I 
IF 

AI 


62  K  = 
AR  (K  ■ 
AI(K, 
K,  JJ) 
K,  JJ) 
K ,  N  )  = 
K ,  N )  = 
1 


1,N 

,K) 
tK) 
-AR(N,K) 

=  AI  (,\|,K) 

7  7 

72 

ltN 

J  J  3 

'J  J) 

=AR{ K,N) 

=  A I  {  K  ,  N ) 

72 


66  1  = 
66  J= 
It  J)  = 

I  V   J  )  = 

600  I 
600  J 
I ,  J  )  = 
I,  J)  = 
702  I 
702  J 
I ,  J  )  = 
It  J)  = 
(  I- J) 
I,  l)  = 
Itlh 
TINUE 


I  3  N 

1  ,N 
AR(  Is 
AIM, 
=  lfM 
=  I«N 
ARff, 

Aid, 

=  lt-M 

=  i;h 

CRU, 

ZVA: 


J)-XR(I)*AI (N+1,J)-XI(I )»AR(N+1 Jj) 


J  ; 


J) 


J) 
J) 

701 .702 
I  )-ALR 
I) -ALI 


1  1 


ASSIGN  911  ''TO  I A 
GO  TO  535 


530 


TO 


A 


IB 
IC 


914  ISL=-T 

915  DO  703  1  =  },;-: 
I I)=0. 

XI ( i )=0. 
ASSIGN  753  70 
ASSIGN  7  04  70 
GO  70  530 
ASSIGN  525  TO  IC 
CO  TO  92 
DO  920  I  =  1,m 
T  TWO  INSTRUCTIONS  FROM 
EIGENVECTORS  XR  AND  XI 
GRi  I , LOO?)  =  XR(  I) 
GI ( I, LOOP)  =  XI!  T  } 
ASSIGN  40  TO  IB 
LOOP  =  LOOP  -f-  1 
IF(N-l)  921,67,523 
ALR  =  ARM,1> °',D^ 
ALI=AI(1 ,1) 
VALUES! 1 -LOOP)  =  ALR 


703 


916 

704 

NE 


920 


525 
67 


TRANSFORMATION  MATRICES  GR  AND  GI  FROM 


VALUESC2,LC0P)  = 
N=0 

GO  TO  700 
GO  TO  4000 


ALI 


3 

1CC0 

1 


3000 
4000 

990 
991 
992 


2 

10 


11 
12 
13 

1U 


15 

20 
30 


32 

33 

35 
36 
34 
1*0 


41 

42 

43 


5  0 
75 


f'l  A  — 

SAVE 
VALU 

VALU 

SAVE 

GfUK 
GR(K 
E 
GKK 
G I  ( K 
CONT 
GO  T 

FORM 

SUBR 

DIME 
OG  1 
DO  1 


000 

=  V 

I 

E 
MX-J 
ES(  1 
ES(1 

=  V 

ES(2 
000 

,J) 

=  G 


10.  **  20 

-■  VALUE  S(  1,1))  3 
ALUESt 1,1) 


3, 1000 


I 

,M 

»J 
AL 
tH 
J 
K  = 
R( 


I( 


,MX)  = 

,J3  = 
INUE 


1  .4 
X}_ 

UES 
X)  = 
)  = 
1,M 
K  I  M 
GR( 
SAV 
K ,  M 
GI  ( 
SAV 


000,1 

=  VALUES* 1 , J) 

I  AX 
( 2 , MX  ) 

VALUES (2, J) 

SAVE 

X) 

KtJ  ) 
E 

X) 

Kf  JJ 
E 


T  99 

AT(15H  NO  CONVERGENCE) 


1 


1  X( 


DO  2 
X  ( K , 

DO  3 
K?  =  0 
Z-0. 
DO  1 

ir-n 

K  P  =  K 
CONT 

IF(L 
CO  1 
Z  =  A[ 


A( 


t~ , 


OUT  I 

N  S  1 0 

1  =  1 

J=l 

J)=0 

;< }  =  i 

4  L- 

0 

2  K= 
-ABS 
SF(A 

[NUE 
-KP) 

4  J  = 
L,J) 
J)=A 
,J)  = 

5  J= 
L  ,  J ) 
J)~X 
,  J  )  = 
BSF{ 
-N)3 
L  +  l 

6  K  = 
(K;l 

3  J  = 
J)=A 

5  J  = 

J)=X 

3  1  = 
+  1-1 
3  J  = 
0 

-N) 

=  11  + 

-2  K= 

A{  II 

J)  = 
1 

0  75 
KER--2 

CONTINUE 

END 

END 


DO  1 
Z  =  X( 

ML, 

IF  (A 
IF(L 
LP1  = 

DO  3 
IF  (A 
R  A  T I 

no  3 

A(K, 

X(K, 
CCNT 
CONT 

1 1-;-; 
DO  4 
S=^0U 
IFl  I 
IIP! 

s=s+, 

XI  II 
KER; 

GO 


NE    GAUSS3(N,E?,A,X,K 
N  A (40, 40),  X( 40,40) 

,N- 

•  Q 

•  N  ■ 
*Q 

1 9.W 


L,N 
F(A(K,L)))11,12,12 

Uy-L)  ) 

13^0,20 

(KP, J) 
7 

1 » N 

(KP, J) 

Z 

A(L,L)  )-EP) 50,50,30 

1,34,34 

LP!  ,N 

) )32,36,32 

K,L)/A(L,L) 

L  P 1  ,  N 

(K, J)-RATIO*A(L, J) 

ltN 

IK, J)-RATIO»X(L, J) 


1,N 

1  i  N 

41,43,43 

1 1 P  1  ,  N 

.  K )  *  X  ( K  ,  J ) 

( XC 1 1 ,J)-S)/A( II ,11) 


ER) 


430 


APPENDIX  II 

GRAPHS 
DF  &ZW11ATED  SYSTEM  TRAJECTORIES 


44 


47 
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